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ABSTRACT 

We performed high-resolution simulations of a sample of 14 galaxy clusters that span a mass 
range from 5 x 10 13 /z~' M to 2 x 10 15 /z 1 M to study the effects of cosmic rays (CRs) on 
thermal cluster observables such as X-ray emission and the Sunyaev-ZeFdovich effect. We 
analyse the CR effects on the intra-cluster medium while simultaneously taking into account 
the cluster's dynamical state as well as the mass of the cluster. The modelling of the cosmic ray 
physics includes adiabatic CR transport processes, injection by supernovae and cosmological 
structure formation shocks, as well as CR thermalization by Coulomb interaction and catas- 
trophic losses by hadronic interactions. While the relative pressure contained in CRs within 
the virial radius is of the order of 2 per cent in our non-radiative simulations, their contribu- 
tion rises to 32 per cent in our simulations with dissipative gas physics including radiative 
cooling, star formation, and supernova feedback. The relative CR pressure rises towards the 
outer regions due to a combination of the following effects: CR acceleration is more efficient 
at the peripheral strong accretion shocks compared to weak central flow shocks, adiabatic 
compression of a composite of CRs and thermal gas disfavours the CR pressure relative to 
the thermal pressure due to the softer equation of state of CRs, and CR loss processes are 
more important at the dense centres. Interestingly, in the radiative simulations the relative CR 
pressure reaches high values of the order of equipartition with the thermal gas in each cluster 
galaxy due to the fast thermal cooling of gas which diminishes the thermal pressure support 
relative to that in CRs. This also leads to a lower effective adiabatic index of the composite gas 
that increases the compressibility of the intra-cluster medium. This effect slightly increases 
the central density, thermal pressure and the gas fraction. While the X-ray luminosity in low 
mass cool core clusters is boosted by up to 40 per cent, the integrated Sunyaev-Zel'dovich 
effect appears to be remarkably robust and the total flux decrement only slightly reduced by 
typically 2 per cent. The resolved Sunyaev-Zel'dovich maps, however, show a larger variation 
with an increased central flux decrement. 

Key words: radiation mechanisms: general - cosmic rays - galaxies: cluster: general - cool- 
ing flows - large-scale structure of Universe - X-rays: galaxies: clusters 



1 INTRODUCTION 

Cosmic ray protons (hereafter referred to as CRs) play a decisive 
role within the interstellar medium of our own Galaxy. CRs and 
magnetic fields each contribute roughly as much energy and pres- 
sure to the galactic ISM as the thermal gas does, CRs trace past 
energetic events such as supernovae, and they reveal the underlying 
structure of the baryonic matter distribution through their interac- 
tions. In contrast, it is unknown how much pressure support is pro- 
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vided by CRs to the thermal plasma of clusters of galaxies. A sub- 
stantial CR pressure contribution might have a major impact on the 
properties of the intra-cluster medium (ICM) and potentially mod- 
ify thermal cluster observables such as the X-ray emission and the 
Sunyaev-Zel'dovich (SZ) effect (Sunyaev & Zeldovich 1972). This 
effect arises since photons from the cosmic microwave background 
experience inverse Compton up-scattering by thermal electrons of 
the ICM. CR induced modifications might imprint themselves in 
changes of cluster scaling relations or modify their intrinsic scat- 
ter and thus change the effective mass threshold of X-ray or SZ 
surveys. Neglecting such a CR component in reference simulations 
might introduce biases in the determination of cosmological param- 
eters. Finally, simulating realistic CR distributions within galaxy 
clusters enables us to predict the non-thermal cluster emission re- 
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suiting from CR interactions with thermal gas protons. The result- 
ing model signatures of the expected radio synchrotron and y-ray 
emission can then be tested with the upcoming new generation of 
y-ray instruments (e.g. imaging atmospheric Cerenkov telescopes 
and the GLAST 1 satellite) and radio telescopes (e.g. LOFAR 2 and 
extended Very Large Array 3 ). 

A thorough dynamical treatment of CR transport and loss pro- 
cesses in numerical simulations of cosmological large-scale struc- 
ture formation has not yet been attempted due to the very complex 
CR physics involved, despite the potential dynamical importance of 
CRs. There have been pioneering efforts to study the CR acceler- 
ation into the intergalactic medium (IGM) on large scales (Miniati 
et al. 2001; Kang et al. 2002; Ryu & Kang 2003, 2004; Kang & 
Jones 2005). However, these codes were quite limited in their adap- 
tive resolution capability and neglected dissipative gas physics in- 
cluding radiative cooling, star formation, and supernova feedback. 
There have also been numerical modelling of discretised CR en- 
ergy spectra on top of grid-based cosmological simulations (Miniati 
2001), but so far these implementations neglected the hydrody- 
namic pressure of the CR component and were unable to resolve 
the observationally accessible, dense central regions of clusters. 

First important results underscoring the importance of CR 
populations in clusters were obtained in non-radiative cosmolog- 
ical simulations and predicted relative CR pressure contributions 
°f Pcs.lP\h ~ 20 per cent up to 50 per cent (Miniati et al. 2001), 
which also motivated suggestions that CRs may slow-down clus- 
ter cooling flows (Cen 2005). Since the CR acceleration depends 
crucially on the shock strength of cosmological structure forma- 
tion shocks, there has been a series of efforts undertaken to develop 
schemes that are capable of measuring the Mach numbers in cos- 
mological simulations (Miniati et al. 2000; Ryu et al. 2003; Pfrom- 
mer et al. 2006). To allow studies of the dynamical effects of CRs 
in radiatively cooling galactic and cluster environments, we have 
developed a CR formalism that is based on smoothed particle hy- 
drodynamical representation of the equations of motion. We intro- 
duced a number of approximations to reduce the complexity of the 
problem while still being able to capture as many physical proper- 
ties and peculiarities of CR populations as possible. In our model, 
the emphasis is given to the dynamical impact of CRs on hydrody- 
namics, and not on an accurate spectral representation of the CRs. 
The guiding principles are energy, pressure, and particle number 
conservations, as well as adiabatic invariants. Non-adiabatic pro- 
cesses are mapped onto modifications of these principles (EnBlin 
et al. 2006; Jubelgas et al. 2006; Pfrommer et al. 2006). 

A simplified overview over the relevant physical processes in 
galaxy clusters is given in Fig. 1. The upper central part shows 
standard processes which are usually considered in simulations. 
Radiative cooling of the gas leads eventually to star formation in the 
densest regions that exceed a certain density threshold. This hap- 
pens in the central cluster regions and within substructures leading 
to individual galaxies. Once the nuclear energy has been used up, 
massive stars explode in supernovae that drive strong shock waves 
into the ambient interstellar medium (ISM) which resupply thermal 
and turbulent energy. On larger scales, structure formation shock 
waves dissipate gravitational energy associated with hierarchical 
clustering into thermal energy of the gas, thus supplying the intra- 
cluster medium (ICM) with entropy and thermal pressure support. 
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Cluster observables: Physical processes in clusters: 
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Figure 1. Overview over the relevant physical processes in galaxy clus- 
ters. The right-hand side shows the interplay of different physical processes 
while the left-hand side shows observables that inform about the properties 
of cluster and their dynamical state. Gain processes are denoted in green 
and loss respectively redistribution processes are denoted in red. 



There are three main observables associated with these processes: 
the hot ICM emits thermal bremsstrahlung radiation with an X- 
ray luminosity that depends on the square of the electron number 
density. The amplitude of the Sunyaev-Zel'dovich effect (Sunyaev 
& Zeldovich 1972) depends on the pressure of the thermal electron 
population integrated along the line-of-sight through the cluster. Fi- 
nally, galaxy spectra probe directly the stellar populations of intra- 
cluster galaxies and indirectly the cluster's potential through their 
velocity dispersion (for reviews see Sarazin 1988; Voit 2005). 

The lower part of Fig. 1 sketches the cosmic ray physics 
within clusters. CRs behave differently compared to the thermal 
gas. Their equation of state is softer, they are able to travel actively 
over macroscopic distances, and their energy loss time-scales are 
typically larger than the thermal ones. Besides thermalization, col- 
lisionless shocks are also able to accelerate ions of the high-energy 
tail of the Maxwellian through diffusive shock acceleration (for re- 
views see Drury 1983; Blandford & Eichler 1987; Malkov & O'C 
Drury 2001). These energetic ions are reflected at magnetic irregu- 
larities through magnetic resonances between the gyro-motion and 
waves in the magnetised plasma and are able to gain energy in mov- 
ing back and forth through the shock front. This acceleration pro- 
cess typically yields a cosmic ray proton (CR) population with a 
power-law distribution of the particle momenta. CRs are acceler- 
ated on galactic scales through supernova shocks while they are 
injected by structure formation shock waves on even larger scales 
up to tens of Mpc. 

The ubiquitous cosmic magnetic fields couple the otherwise 
dynamically independent ingredients like the ICM plasma, and the 
CR gas into a single, however complex fluid. Magnetic fields pre- 
vent charged relativistic particles to travel macroscopic distances 
with their intrinsic velocity close to the speed of light. Instead, 
the particles gyrate around, and travel slowly along magnetic field 
lines. Occasionally, they get scattered on magnetic irregularities. 
On macroscopic scales, the transport can often be described as a 
diffusion process if the gyro-radius can be regarded to be small. 
Thus, CR diffusion redistributes the CR energy density macroscop- 
ically. Thermal heat conduction is an analogous process that real- 
locates the thermal energy of the ICM. The CR energy reservoir 
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suffers two main loss processes: (1) CR energy is transferred into 
the thermal energy reservoir through individual electron scatterings 
in the Coulomb field of the CR particle as well as by small momen- 
tum transfers through excitations of quantised plasma oscillations. 
We refer to the sum of both effects as Coulomb losses. (2) Provided 
the CR momentum exceeds the threshold p =s 0.8 GeV/c for the 
hadronic reaction with ambient protons, they produce pions which 
decay into secondary electrons (and neutrinos) and y-rays: 

71* -> pr + v^/v^, -> e ± + v e lv e + + 
n° -» 2y. 

Only CR protons above this kinematic threshold are therefore vis- 
ible through their decay products via radiative processes, making 
them directly observationally detectable. These hadronically pro- 
duced relativistic electrons and positrons can emit a halo of radio 
synchrotron emission in the presence of ubiquitous intra-cluster 
magnetic fields. However, there are other processes that acceler- 
ate relativistic electrons such as resonant and non-resonant interac- 
tions of mildly relativistic electrons with magneto-hydrodynamic 
waves or turbulent spectra. Since the distribution of magnetic field 
strengths with cluster radius is also not well known, radio syn- 
chrotron emission alone has limited predictive power. Future y-ray 
satellites should be able to detect the associated hadronically in- 
duced y-ray emission resulting from neutral pion decay and allow 
stronger conclusions on the parent CR population in clusters. 

So far, we have neglected feedback from active galactic nuclei 
(AGN) in our simulations despite its importance (for first models of 
the thermal feedback within cosmological simulations, see Sijacki 
& Springel 2006). Gravitational energy associated with the accre- 
tion onto super-massive black holes is converted into large-scale 
jets and eventually dissipated into thermal and CR energy. A reli- 
able model would be required to self-regulate and respond to the 
instantaneous cooling rate. 

This paper studies the interplay of thermal gas and CRs and 
their effect on the observables of the thermal gas such as X-ray 
emission and the Sunyaev-Zel'dovich effect. Two companion pa- 
pers are studying directly the CR related multi-frequency observ- 
ables and the systematic bias and scatter of cluster scaling relations 
that are associated with CRs and how these modifications effect 
the determination of cosmological parameters. The outline of the 
paper is as follows. Section 2 describes the general setup of the 
simulations, the different physical processes which we simulated, 
and our cluster sample. In Sect. 3, we present the results on the 
relative CR pressure integrated over the clusters, and show the pro- 
files of thermodynamical and CR-related quantities separately for 
non-radiative and radiative simulations. Both types of simulation 
models will be concluded by an emerging picture of CR-induced 
changes to the ICM. Section 3.4 analyses physical processes related 
to CRs in detail and complements the previous section. Finally, 
Sect. 3.5 studies the consequences of these processes for thermal 
observables such as X-ray emission and the Sunyaev-Zel'dovich 
effect. 



2 THE SIMULATIONS 
2.1 General setup 

Simulations were performed using the "concordance" cosmologi- 
cal cold dark matter model with a cosmological constant (ACDM). 
The cosmological parameters of our model are: Q m = Q DM + fl b = 
0.3, n b = 0.039, fi A = 0.7, h = 0.7, n = 1, and <r s = 0.9. Here, Cl m 



denotes the total matter density in units of the critical density for ge- 
ometrical closure today, p crit = 3Hq/(StzG). fit, and Q A denote the 
densities of baryons and the cosmological constant at the present 
day. The Hubble constant at the present day is parametrized as 
Hq= 100/j km s~ Mpc~ , while n denotes the spectral index of the 
primordial power-spectrum, and <x 8 is the rms linear mass fluctua- 
tion within a sphere of radius 8 /r'Mpc extrapolated to z = 0. This 
model yields a reasonable fit to current cosmological constraints. 

Our simulations were carried out with an updated and ex- 
tended version of the distributed-memory parallel TreeSPH code 
GADGET-2 (Springel 2005; Springel et al. 2001) that includes self- 
consistent cosmic ray physics (EnBlin et al. 2006; Jubelgas et al. 
2006; Pfrommer et al. 2006). Gravitational forces were computed 
using a combination of particle-mesh and tree algorithms. Hydro- 
dynamic forces are computed with a variant of the smoothed parti- 
cle hydrodynamics (SPH) algorithm that conserves energy and en- 
tropy where appropriate, i.e. outside of shocked regions (Springel 
& Hernquist 2002). 

We have performed high-resolution hydrodynamic simula- 
tions of the formation of 14 galaxy clusters. The clusters span a 
mass range from 5 x 10 13 h' 1 M Q to 2 x 10' 5 h' 1 M Q and have orig- 
inally been selected from a low-resolution dark-matter- (DM)-only 
simulation (Yoshida et al. 2001) with box-size 479 /T 1 Mpc of a 
flat ACDM model with the parameters mentioned above. Using the 
'zoomed initial conditions' technique (Katz & White 1993; Tormen 
et al. 1997), we then re-simulated the clusters with higher mass and 
force resolution by adding short-wavelength modes within the La- 
grangian regions in the initial conditions that will evolve later-on 
into the structures of interest (using initial conditions from Dolag 
et al. 2005). In these high-resolution regions, dark matter particles 
of the parent simulation are split into a dark matter and gaseous 
part, with the mass ratio reflecting the value of the cosmic baryon 
fraction. Force and mass resolution are then gradually degraded in 
regions larger than at least 3-5 virial radii from the clusters to limit 
the computational cost while providing a correct representation of 
the large scale tidal gravitational field. Using an iterative method, 
the size of the cluster region that is re-simulated at high-resolution 
was always chosen to be large enough to prevent any contamination 
with heavy dark matter particles. 

In high-resolution regions, the dark matter particles had 
masses of m DM = 1.13 x 10 9 /r'M o and SPH particles m gas = 
1.7 x 10 8 h' 1 M Q so each individual cluster is resolved by 8 x 10 4 to 
4 x 10 6 particles, depending on its final mass. For details of cluster 
properties, see Sect. 2.3. The SPH densities were computed from 48 
neighbours, allowing the SPH smoothing length to drop at most to 
half of the value of the gravitational softening length of the gas par- 
ticles. This choice of the SPH smoothing length leads to our mini- 
mum gas resolution of approximately 8 x 10 9 hr l M Q . For the initial 
redshift we chose 1 + zm = 60. The gravitational force softening 
was of a spline form (e.g., Hernquist & Katz 1989) with a Plum- 
mer equivalent softening length that is assumed to have a constant 
comoving scale down to z = 5, and a constant value of 5 /r'kpc in 
physical units at later epochs. 

2.2 The models 

For each galaxy cluster we ran six different simulations separated 
into two sets of three simulations each (cf. Table 1). One set uses 
only non-radiative gas physics while the other set consists of three 
simulations that include radiative cooling and star formation. The 
first simulation of the non-radiative set is our reference simulation 
with non-radiative physics only, i.e. the gas is transported adiabat- 
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Table 1. : Different physical processes used in our simulation models: 



Simulated physics' ' non-radiative simulation models' radiative simulation models' ': 

reference mod. shock-CR mod. simple shock-CR mod. reference mod. shock-CR mod. complete CR mod. 



thermal shock heating •/ 


/ 




/ 


/ 


/ 


radiative cooling 






/ 


/ 


/ 


star formation 






/ 


/ 


/ 


Coulomb CR losses 


■/ 


/ 




/ 


/ 


hadronic CR losses 


/ 


/ 




/ 


/ 


shock-CRs with f = const' 2 ' 




/ 








shock-CRs with <T(M)' 2 > 


/ 






/ 


/ 


supernova-CRs 










/ 



Notes: 

(1) This table serves as an overview over our simulated models. The first column shows the simulated physics, the next three columns specify our non-radiative 
simulation models, while the last three columns show our radiative simulation models. 

(2) We employ two different CR injection schemes at structure formation shocks: The realistic one uses a Mach number dependent energy injection efficiency 
((M) while the simplified scheme uses a constant energy injection efficiency, ( = const. The latter scheme is a simplification that exaggerates the CR effects 
for better visibility. For further details and references please refer to Sec. 2.2. 



ically unless it experiences structure formation shock waves that 
supply the gas with entropy and thermal pressure support. The sec- 
ond simulation of the set allows for a simplified model of CR accel- 
eration at structure formation shocks with a constant energy injec- 
tion efficiency of £ = £cR,inj/£diss = 0.5, independent of the shock 
strength, while reducing at the same time the amount of thermal 
dissipation to ensure energy conservation at the shock. Here, ScRjnj 
denotes the injected CR energy density and e diss the total dissipated 
energy at the shock. This model represents an extreme simplifica- 
tion that exaggerates the CR effects for better visibility. All CR 
models include Coulomb and hadronic loss processes. The third 
simulation of the non-radiative set uses a more realistic model for 
the CR acceleration at structure formation shocks which is based 
on the thermal leakage model (e.g., Ellison & Eichler 1984; Kang 
& Jones 1995; EnBlin et al. 2006). It employs the fact that the 
most energetic ions in the exponential tail of the Maxwellian of 
the post-shock regime are able to participate in the process of diffu- 
sive shock acceleration and injects a power-law in momentum that 
depends on the instantaneous shock strength in the simulation (for 
details, cf. Pfrommer et al. 2006). With this thermal leakage model, 
we are able to derive a CR energy injection efficiency, £(M), for the 
diffusive shock acceleration process. It depends only on the Mach 
number of the shock, M = v s / c i (shock velocity v s divided by 
pre-shock sound speed ci), the post-shock temperature T2, and a 
saturation value of the CR energy injection efficiency at high val- 
ues of the Mach number, f max = 0.5 as shown in Fig. 2 (Ryu et al. 
2003; EnBlin et al. 2006). 

The first simulation of our radiative set is a reference simula- 
tion with radiative cooling and star formation, but no CR physics. 
The second radiative simulation accounts for CR acceleration at 
structure formation shocks and allows for all CR loss processes. Fi- 
nally, the last simulation additionally assumes that a constant frac- 
tion £sn = 0.3 of the kinetic energy of a supernova ends up in the 
CR population, which is motivated by TeV y-ray observations of a 
supernova remnant that find an energy fraction of £ SN ^ 0.1 - 0.3 
when extrapolating the CR distribution function (Aharonian & et 
al. 2006). We choose a maximum value for this supernova energy 
efficiency owing to the large uncertainties and our aim to bracket 
the realistic case with the two radiative CR simulations. 

Radiative cooling was computed assuming an optically thin 
gas of primordial composition (mass-fraction of X = 0.76 for hy- 
drogen and 1 - X = 0.24 for helium) in collisional ionisation 



Mach number M 




Figure 2. CR energy injection efficiency for the diffusive shock accelera- 
tion process. Shown is the CR energy injection efficiency £ = £cR,inj/£diss 
for the three post-shock temperatures i^/keV = 0.01, 0.3, and 10. For this 
illustrative figure, we consider only CRs above the kinematic momentum 
threshold threshold = Pcr/("i p c) = 0.83 of the hadronic CR p-p interac- 
tion which are able to produce pions that successively decay into secondary 
electrons, neutrinos, and y-rays. We choose equipartition between injected 
CR energy and dissipated thermal non-relativistic energy as plausible satu- 
ration value of the CR energy injection efficiency, i.e. £ mllx = 0.5 (Ryu et al. 
2003; EnBlin et al. 2006). 



equilibrium, following Katz et al. (1996). We have also included 
heating by a photo-ionising, time-dependent, uniform ultraviolet 
(UV) background expected from a population of quasars (Haardt 
& Madau 1996), which reionizes the Universe at z - 6. The ef- 
fect of a photo-ionising background is that of preventing the col- 
lapse of overdense regions of gas and inhibiting subsequent star 
formation within halos of sub-L* galaxies (e.g., Benson et al. 2002), 
thus having a secondary effect on the resolution of our simulation. 
If not counteracted by any feedback process, radiative cooling in 
high-resolution simulations is well known to give rise to a central 
cooling flow that overproduces the amount of stars not only in the 
group/cluster overdense environment but also in the average envi- 
ronment (e.g., Balogh et al. 2001; Borgani et al. 2004, 2006, and 
references therein). As discussed by Springel & Hernquist (2003), 
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their multiphase ISM model with supernovae feedback alone does 
however not fully resolve this problem, despite its ability to reg- 
ulate the consumption of cold gas into stars within the ISM. The 
reason for this is because the cooling rates within halos remain es- 
sentially unaffected in the model, i.e. the supply of gas to the dense 
star-forming ISM is largely unchanged, while by construction the 
phases of the ISM remain coupled to each other, preventing baryons 
to leave the ISM (except for dynamical effects like gas stripping in 
mergers). Being mainly interested in global properties of the ICM, 
we do not include the metallicity effect in the cooling function due 
both to the approximate treatment of metal generation and diffusion 
and in order not to exaggerate the cooling problem of the cluster 
cores since we also do not follow the AGN feedback that seems to 
be essential to reliably model the dense cluster centres (e.g., Sijacki 
& Springel 2006). 

Star formation is treated using the hybrid multiphase model 
for the interstellar medium introduced by Springel & Hernquist 
(2003). In short, the ISM is pictured as a two-phase fluid consisting 
of cold clouds that are embedded at pressure equilibrium in an am- 
bient hot medium. The clouds form from the cooling of high den- 
sity gas, and represent the reservoir of baryons available for star 
formation. When stars form, according to a Salpeter initial mass 
function (IMF) (Salpeter 1955), the energy released by supernovae 
heats the ambient hot phase of the ISM, and in addition, clouds in 
supernova remnants are evaporated. These effects establish a tightly 
self-regulated regime for star formation in the ISM. In practice, 
the scheme is numerically implemented as a sub-resolution model, 
i.e. cold clouds are not resolved individually. Instead, only their to- 
tal mass fraction in each element of the ISM is computed, otherwise 
they are treated in a stochastic fashion with their collective effect 
on the ISM dynamics being described by an effective equation of 
state. The numerical implementation of this multiphase model de- 
scribes each gas particle as composed by a hot component, having 
its own mass and density, and a cold neutral component. The rela- 
tive amount of these two phases is determined by the local value of 
gas density and temperature. 

Cosmic ray physics was computed by using a new formula- 
tion that follows the most important injection and loss processes 
self-consistently while accounting for the CR pressure in the equa- 
tions of motion (EnBlin et al. 2006; lubelgas et al. 2006; Pfrom- 
mer et al. 2006). We refer to these papers for a detailed descrip- 
tion of the formalism, providing here only a short summary of the 
model. In our methodology, the non-thermal cosmic ray population 
of each gaseous fluid element is approximated by a simple power 
law spectrum in particle momentum, characterised by an ampli- 
tude, a low-momentum cut-off, and a fixed slope a = 2.3. Adia- 
batic CR transport processes such as compression and rarefaction, 
and a number of physical source and sink terms which modify the 
cosmic ray pressure of each particle are modelled. The most impor- 
tant sources considered are injection by supernovae (in our radia- 
tive simulations) and diffusive shock acceleration at cosmological 
structure formation shocks, while the primary sinks are thermaliza- 
tion by Coulomb interactions, and catastrophic losses by hadronic 
interactions. 

2.3 Simulated clusters 

We chose our sample of 14 galaxy clusters such that they span a 
wide range in mass (5 x 10 13 /r' M Q to 2 x 10 15 /r' M Q ) and show 
a variety of dynamical states ranging from relaxed cool core (CC) 
clusters to violent merging clusters (cf. Table 2). Since the CR in- 
jection crucially depends on the strength of shock waves we expect 



Table 2. : Cluster sample 



Cluster 


sim.'s 


dyn. stated 


w<2) 
200 
[h~ l Mq] 


n(2) 
200 

[/r 1 Mpc] 


kT Q) 
200 
[keV] 


1 


g8a 


CC 


1.8 x 10 15 


2.0 


13.1 


2 


gla 


CC 


1.3 x 10 15 


1.8 


10.6 


3 


g72a 


PostM 


1.1 x 10 15 


1.7 


9.4 


4 


g51 


CC 


1.1 x 10 15 


1.7 


9.4 


5 


gib 


M 


3.7 x 10 14 


1.2 


4.7 


6 


g72b 


M 


1.5 x 10 14 


0.87 


2.4 


7 

1 


glc 


1V1 


1 A n/ 1 ni4 


U.04 


Z.j 


8 


g8b 


M 


l.Ox 10 14 


0.76 


1.9 


9 


gld 


M 


9.2 x 10 13 


0.73 


1.7 


10 


g676 


CC 


8.8 x 10 13 


0.72 


1.7 


11 


g914 


CC 


8.5 X 10 13 


0.71 


1.6 


12 


gle 


M 


6.4 X 10 13 


0.65 


1.3 


13 


g8c 


M 


5.9 X 10 13 


0.63 


1.3 


14 


g8d 


PreM 


5.4 X 10 13 


0.61 


1.2 



Notes: 

(1) The dynamical state has been classified through a combined criterion 
invoking a merger tree study and the visual inspection of the X-ray bright- 
ness maps. The labels for the clusters are M-merger, PostM-post merger 
(slightly elongated X-ray contours, weak cool core region developing), 
PreM-pre-merger (sub-cluster already within the virial radius), CC-cool 
core cluster with extended cooling region (smooth X-ray profile). 

(2) The virial mass and radius are related by M&(z) = \n Ap a n(z)R A , 
where A = 200 denotes a multiple of the critical overdensity p cr it(z) = 
3H(z) 2 /(SnG). 

(3) The virial temperature is defined by kT & = GM A fim p /(2R A ), where p 
denotes the mean molecular weight. 

the impact of CR effects to depend on both, the cluster size and its 
dynamical state. To investigate this dependence, we re-simulated 
three isolated clusters (cluster 4, 10, and 1 1) and three super-cluster 
regions which are each dominated by a large cluster (cluster 1, 2, 
and 3) and surrounded by satellite clusters (cluster 5 to 9 and 12 to 
14). 

We analysed the clusters with a halo-finder based on spherical 
overdensity followed by a merger tree analysis in order to get the 
mass accretion history of the main progenitor. We also produced 
projections of the X-ray emissivity at redshift z = in order to 
get a visual impression of the cluster morphology. The dynamical 
state of a cluster is defined by a combined criterion: (i) if the clus- 
ter did not experience a major merger with a progenitor mass ratio 
1:3 or larger after z = 0.8 (corresponding to a look-back time of 
~ 5 h~ l Gyr) and (ii) if the visual impression of the cluster's X- 
ray morphology is relaxed, it was defined to be a cool core cluster. 
The merging clusters are subdivided into currently merging clus- 
ters, pre-mergers that show a merging sub-cluster already within 
the virial radius, and post-mergers with slightly elongated X-ray 
contours in combination with a weak cool core region developing. 
While this classification may not always be unambiguous in indi- 
vidual cases, it provides robust samples of clusters in these differ- 
ent dynamical states. Interestingly, we note that each satellite clus- 
ter within our three super-cluster regions shows signs of merging 
activity, presumably triggered by the high ambient density that is 
necessary for the development of a super-cluster region. 

The spherical overdensity definition of the virial mass of the 
cluster is given by the material lying within a sphere centred on a lo- 
cal density maximum, whose radial extend R A is defined by the en- 
closed threshold density condition M(< R a )I{AjiR?J2>) = p tbre& . We 
chose the threshold density p t h rcs (z) = Ap cr (z) to be a multiple A = 
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Table 3. : Volume averaged values for (Xcr) = (Pcs)l(Pth) AND EcR/Eth within the virial radius: 



relative CR quantity 


sample' 1 ' 


radiative, 


radiative, 


non-radiative, 


non-radiative, 






shock-CRs <2) 


shock- & SNe-CRs (2) 


shock-CRs (2) 


shock-CRs with f = const' 2 ' 


(Xcr) = <P C R>/<P<h> 


small CC 


0.045 


0.216 


0.022 


0.321 




large CC 


0.034 


0.068 


0.015 


0.369 




merger 


0.161 


0.414 


0.022 


0.249 




all 


0.120 


0.321 


0.021 


0.272 


Ecr/Eo, 


small CC 


0.083 


0.358 


0.037 


0.541 




large CC 


0.057 


0.114 


0.025 


0.609 




merger 


0.266 


0.684 


0.036 


0.410 




all 


0.199 


0.531 


0.034 


0.449 



Notes: 

(1) The different cluster samples consist of small cool core (CC) clusters (cluster 10, 11), large cool core clusters (cluster 1, 2, 4), merging clusters (other 
clusters), and all 14 clusters. 

(2) The third column derives from our radiative simulations that account only for structure formation CRs while the values of the forth column correspond to 
the radiative simulations that take additionally supernova CRs into account. The values in the fifth column derive from non-radiative simulations that accelerate 
CRs at structure formation shocks using a Mach number dependent energy injection efficiency ((M), while the last column derives from the scheme with a 
constant energy injection efficiency, f = const. 



200 of the critical density of the universe p cr (z) = 3H(z) 2 /(SttG). 
For a better comparison with observations and reasons laid out 
in detail by Evrard (2004), we assume a constant A = 200 al- 
though some treatments employ a time-varying A in cosmologies 
with fi m + 1 (Eke et al. 1996). In the reminder of the paper, we use 
the terminology R vk instead of R 2 oo- 



3 RESULTS 

3.1 Relative CR pressure 

Which galaxy clusters are mostly affected by CRs? And does cool- 
ing and star formation change the generic picture of CRs in clus- 
ters? To answer these two questions it is instructive to consider the 
CR pressure P C r relative to the thermal pressure P th integrated over 
the cluster. This is a good measure of the importance of dynamical 
effects of CRs on the ICM. Table 3 shows the volume averaged val- 
ues for (Xcr) = (PcR)/(Pth) within the virial radius R v „ and the 
ratio of the CR-to-thermal energy within the galaxy cluster. 

The main findings are the following. (1) Among the non- 
radiative simulations, (X C r) is much smaller ((X CR ) ~ 2 per cent) 
in our improved scheme where CRs are accelerated at structure for- 
mation shocks using Mach number dependent energy injection ef- 
ficiency f (M) compared to the simplified scheme where CRs have 
been injected with a high constant energy injection efficiency of 
( = 0.5. This shows the importance of a Mach number depen- 
dent CR injection scheme to study CR effects within the ICM. 
(2) Neglecting for the moment the simplified scheme, we notice 
that (X C r) is much larger in our radiative simulations compared to 
our non-radiative simulations for reasons that will be analysed in 
the next sections. (3) Among our radiative simulations, the com- 
plete model that accounts for CRs from structure formation shocks 
on large scales and supernova shocks within the ISM of individ- 
ual galaxies shows a larger relative CR pressure compared to the 
model that only accounts for structure formation CRs. Note that 
the complete model opens simply another channel for CR injection. 
(4) Comparing different dynamical states of clusters, the relative 
CR pressure contribution is largest in merging clusters, followed 
by small cool core clusters, and large cool core clusters. Strong 
merger shock waves efficiently inject CRs at the external cluster 



regions and mix the highly CR-enriched intergalactic medium out- 
side clusters with the ICM. An effect that boosts the relative CR 
pressure during a merger. Weak virialisation shocks traversing the 
cluster after the merger thermalize random gas motions and yield a 
decrease of (X C r}- (5) The volume averages for X C r are dominated 
by the outer regions due to the rising profile of X C r with increas- 
ing radius. In other words, the absolute value of (X C r) depends on 
the definition of the virial radius. A larger virial radius Z?ioo > ^200 
would thus lead to an increase of (X C r). 

The low level of (Xcr} in our non-radiative simulations dis- 
agrees with the substantially larger predictions of (Xcr) ~ 20 per 
cent to 50 per cent by Miniati et al. (2001). This could be due to 
a number of reasons. Miniati et al. (2000) identified shocks with 
Mach numbers in the range 4 < M < 5 as the most important in 
thermalizing the plasma. In contrast, Ryu et al. (2003) and Pfrom- 
mer et al. (2006) found that the Mach number distribution peaks in 
the range 1 < M < 3. Since diffusive shock acceleration of CRs 
depends sensitively on the Mach number (cf. Fig. 2), this implies 
a more efficient CR injection in the simulations by Miniati et al. 
(2001). Secondly, CR injection according to the "thermal leakage" 
model in the implementation of Miniati (2001) has no injection 
limiter for the saturated stage of non-linear CR acceleration which 
might have led to an overestimate of the CR injection for strong 
shocks. Finally, the grid-based cosmological simulations have been 
performed in a cosmological box of side-length 50 h' 1 Mpc with 
a spatial resolution of 200 /T 1 kpc, assuming an Einstein-de Sitter 
cosmological model (Miniati et al. 2001). The lack of resolution 
inside the clusters in the grid-based approach does not allow one 
to properly resolve the adiabatic compression of the composite of 
CRs and thermal gas. As a result, on one hand the disfavouring of 
the CR pressure relative to the thermal pressure upon compression 
is not accounted for, and on the other hand, this underestimates CR 
loss processes at the dense cluster centres due the lack of resolution 
there. 

While our simulations have much higher spatial resolution in- 
side clusters and can therefore be expected to be more accurate, one 
has to be aware of resolution limitations in them as well. This ap- 
plies especially to the radiative simulations, which are not able to 
resolve faint cluster galaxies, and are not guaranteed to give already 
a converged value for the total amount of stars forming in the clus- 
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ters. Also, we know that simulated clusters exhibit cooling flows 
that are much stronger than observed in nature, which can amplify 
effects related to the CR compression and slow infall of ICM gas in 
the simulated clusters. However, our main purpose in this paper is 
to study the principle changes in cluster structure brought about by 
the inclusion of CRs, and our resolution is high enough to answer 
this question reliably. 

Taking into account that the two radiative CR simulations are 
chosen to bracket the realistic case, we roughly predict a ratio of the 
CR-to-thermal energy of E C n/E th = 9, 20, and 50 per cent for small 
CC, large CC, and merging clusters. These predictions are in good 
agreement with CR constraints from gamma-ray and radio obser- 
vations (Pfrommer & EnBlin 2003, 2004), considering that these 
authors assumed a constant CR energy density relative to the ther- 
mal gas. 

3.2 Profiles of non-radiative simulations 

Results for the appropriately scaled average profiles of our non- 
radiative simulations are shown in Fig. 3 only for our sample of 
small cool core clusters. The main conclusions drawn from this 
sample are representative for our complete sample, while includ- 
ing the remaining clusters would only increase the cluster to clus- 
ter variance of the profiles, primarily due to the non-sphericity (for 
merging clusters) and the different concentrations (for large cool 
core clusters). The radial coordinate has been scaled with if vir to 
allow comparison between differently sized objects. In all plots of 
the profiles, the error bars represent the standard deviation cr from 
the sample mean of the indicated subsample of clusters. The colour 
scheme for the non-radiative profiles is black for our reference sim- 
ulations without CRs, red for simulations where CRs are injected 
only through shock acceleration using Mach number dependent 
scheme (realistic CR model), and blue for simulations where CRs 
are accelerated at structure formation shocks with a constant injec- 
tion efficiency £ = 0.5 (simplified CR model). One should bare in 
mind that the latter model represents an extreme simplification that 
exaggerates the CR effects for better visibility. 



of a reduced thermal pressure at cluster outskirts and of the in- 
creased central pressure is enhanced in our simplified CR model. 
In the outer regions, the reduced thermal pressure adds up with the 
CR pressure to the total (thermal) pressure in our reference simu- 
lation while the central pressure enhancement reflects the density 
enhancement. The CR pressure shows a flat central plateau and de- 
creases at r > 0. 1-ftyir outwards less steeply compared to the thermal 
pressure. 

This yields a rising profile for the relative CR pressure 
throughout the cluster X CR = Pcr/Poi due to a combination of the 
following effects: CR acceleration is more efficient at the periph- 
eral strong accretion shocks compared to weak central flow shocks, 
adiabatic compression of a composite of CRs and thermal gas dis- 
favours the CR pressure relative to the thermal pressure due to the 
softer equation of state of CRs, and CR loss processes are more 
important at the dense centres. 

3.2.3 Baryon fraction 

In our self-consistent CR model (with varying £(M)), the 
central baryon fraction scaled with the universal value b = 
/bary(< x)£l m /Cl b is increased by ~ 10 per cent and reaches the 
same value at the virial radius of b = 0.9 compared to our refer- 
ence simulation. Note that the baryon fraction is increasing as a 
function of radius which can be traced back to the different nature 
of dark matter and gas. The inner slope of the density profile of col- 
lisionless dark matter does not show any sign of flattening (Navarro 
et al. 2004) while the gas profile is cored due to the isotropic hy- 
drodynamic pressure. This yields to an decreasing baryon fraction 
towards the centre, an effect that has been observed in many works 
(e.g. Ettori et al. 2006, and references therein). 

Remarkably, the baryon fraction in our simplified CR model 
(with £ = 0.5) even exceeds unity at intermediate radii and reaches 
a value at the virial radius of b = 0.95. This is mostly due to the 
enhanced gas compressibility which boosts the density within the 
central regions. 



3.2. 1 Density and temperature profiles 

The gas density has been scaled with fib Petit while the temper- 
ature has been scaled by its virial quantity that is defined by 
kT vl[ = GM vlr p m p /(2R vlr ), where p = 4/(3X H + 1 +4 X H x e ) denotes 
the mean molecular weight, X H = 0.76 is the primordial hydrogen 
mass fraction, and x e is the ratio of electron and hydrogen number 
densities which we dynamically track in our radiative simulations. 
Both the average density and the temperature profile of our realis- 
tic CR model shows only very small differences to our reference 
simulations. The tendency of the enhanced density and the reduced 
temperature in our CR simulations is clearly enhanced in our sim- 
plified CR model. While the density enhancement is largest in the 
centre, the decrease of the cluster temperature is larger at the cluster 
outskirts for reasons that are described below (cf. Sect. 3.2.5). 



3.2.2 Pressure and relative CR pressure profiles 

The lines without symbols denote the thermal pressure, the ones 
with symbols the CR pressure. The CR and thermal pressure have 
been scaled with P vil = kT vir 200 0. b p oAt l{p.m v ). Again, there is 
only a very small difference visible between the thermal pressure 
of our realistic CR model and our reference simulations. The trend 



3.2.4 Effective adiabatic index 

When studying the hydrodynamics of a composite of thermal and 
CR gas, it is appropriate to define an effective microscopic adiabatic 
index by 



Teff 



dlog(P th + /> CR ) 



dlogp 



7th + 7CR ^CR 



(1) 



S 1 + %CR 

where y th = 5/3 is the ideal gas value of the ratio of specific heats 
and the adiabatic exponent of the CR population is defined by 

d log P C r 



Jcr : 



dlogp 



(2) 



while the derivative has to be taken at constant entropy S (EnBlin 
et al. 2006). Depending on the CR momentum distribution func- 
tion, ycR can vary between ycR — * 4/3 (in the ultra-relativistic 
limit) to ycR — > 5/3 (in the non-relativistic limit). One should not 
confuse this effective microscopic adiabatic index of our compos- 
ite fluid that enters the hydrodynamic equations and governs the 
compressibility of the medium with the effective macroscopic adia- 
batic index one can derive from cluster density and pressure profiles 
by means of y = dlog P/(dlog r) [dlogp/(dlogr)]~'. If one com- 
pressed gas with such an effective macroscopic adiabatic index into 
a spherical cluster that has the appropriate density profile (without 
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Figure 3. Average profiles of our sample of small cool core clusters (clusters 10, 1 1) in our non-radiative simulations at redshift z = 0. Here and in the following 
plot, the error bars represent the standard deviation from the sample mean. Shown are the appropriately scaled profiles of the gas density p gas , temperature T, 
thermal and CR pressure (P t h and Pcr, the latter with symbols), relative CR pressure Xqr = Pcn/Pth, baryon fraction b, and the effective adiabatic indices 
y c ff (solid) in comparison with the CR adiabatic indices ycR (dashed). Colour scheme: black: reference simulations without CRs, red: CRs are injected only 
through shock acceleration using our Mach number dependent scheme, blue: CR shock acceleration with constant energy injection efficiency. 
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considering shocks and cooling), one would obtain the simulated 
pressure profiles. 

The effective adiabatic index y eB is shown with a solid line 
while 7cr is shown with a dashed line. As expected, y eff reflects 
the profile of X C r and is only different from the adiabatic index 
7th = 5/3 of the thermal gas in the external cluster regions where 
CRs make up a larger fraction of the pressure support (in partic- 
ular in our simplified model). Interestingly, j C r shows a different 
behaviour with a slightly softer value in the central cluster regions. 
This is due to enhanced Coulomb cooling in denser core region 
which effectively re-thermalizes low-energetic CRs and increases 
the low momentum cutoff of the CR distribution function. This 
leads to a more relativistic CR population with a lower CR pres- 
sure and thus to the observed behaviour of the adiabatic indices. 
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3.2.5 Discussion of CRs profiles in non-radiative simulations 

Structure formation shock waves propagate through the cosmic ten- 
uous gas, which is compressed at the transition layer of the shock 
while a part of the kinetic energy of the incoming gas is dissipated 
into internal energy of the post-shock gas according to the Rankine- 
Hugoniot jump conditions. Considering CR injection at shocks 
through diffusive shock acceleration, a part of the energy that has 
been otherwise completely thermalized is now transferred to the 
CR energy reservoir while employing energy conservation. This ef- 
fect lowers the amount of thermalized energy and thus the temper- 
ature which is observed at the periphery in clusters (cf. Sect. 3.2.1) 
where it is most pronounced due to strong shock waves that take 
place there and inject efficiently CRs (Miniati et al. 2000; Ryu et al. 
2003; Pfrommer et al. 2006). 

The successive build-up of a cluster profile of thermodynamic 
variables (such as entropy or pressure) through mergers and mass 
accretion can be thought of as a combination of two effects. (1) 
Merger shock waves violently increase the cluster's entropy and 
(2) permanent accretion of new shells of mass at the peripheral 
cluster regions exert a pressure that adiabatically compresses the 
central gas. These two effects have important consequences for the 
resulting thermal profiles that are modulated by CR pressure. Since 
CR acceleration is more efficient at the peripheral strong accretion 
shocks compared to weak central flow shocks, this yields a strati- 
fication of the relative CR pressure profile which is increasing out- 
wards. Adiabatic compression of a composite of CRs and thermal 
gas disfavours the CR pressure relative to the thermal pressure due 
to the softer equation of state of CRs. Roughly, one expects the ra- 
tio of X C r = PcR/Pih after adiabatic compression (denoted with a 
prime) to that before compression to scale as 



^CR 



JeL\ 



0.05 



(3) 



for a compression ratio of p'/p ~ 10 4 . Additionally, CR loss pro- 
cesses are more important at the dense cluster centres. Thus, the 
central density has to increase dynamically in our CR simulation 
compared to our reference simulations in order to balance the grav- 
itational pressure and make up for the initially 'missing' thermal 
pressure. Shock heating and adiabatic compression heats the ther- 
mal gas in the cluster centre with a higher rate compared to the ex- 
ternal regions. Thus, the thermal pressure P th = n th kT is increased 
in the cluster centre in our simulations with CRs compared to our 
reference simulations without CRs. 




Figure 4. The top panel shows the cooling times due to Coulomb losses (ris- 
ing solid line) and hadronic dissipation (nearly horizontal line) as a function 
of the spectral cut-off q = p m in/(«ipC). The dot-dashed line gives the total 
cooling time, while the vertical dashed lines marks the asymptotic equilib- 
rium cut-off reached by the CR spectrum when no sources are present. The 
bottom panel shows the cooling time of ordinary thermal gas due to radia- 
tive cooling (for primordial metallicity), as a function of temperature. The 
horizontal dashed line marks the cooling time of CRs with a high momen- 
tum cut-off (q » 1), for comparison. In both panels, the times have been 
computed for a number density of n = 0.01 cirT 3 , a typical number density 
in the centres of galaxy clusters. Note however that the cooling times all 
scale as r oc l/n, i.e. for different densities, only the vertical scale would 
change but the relative position of the lines would remain unaltered. 



3.3 Profiles of radiative simulations 

The larger CR pressure support in our radiative simulations makes 
CR related phenomena more pronounced compared to the non- 
radiative simulations. This crucially depends on the small ratio 
of the thermal to the comparatively long CR cooling time in re- 
gions where the temperature drops below T < 10 7 K. Here, the 
CR cooling time of an aged CR population with an equipartition 
value for the lower momentum cutoff q ~ 1 remains almost con- 
stant, but it decreases for the thermal gas as it is cooling (cf. Fig. 4) 
and becomes progressively shorter than the CR cooling time. This 
equipartition value for q results from a balance of Coulomb cool- 
ing that increases q and leaves the normalisation of the CR dis- 
tribution function almost unchanged and CR cooling by hadronic 
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interactions which decreases the lower momentum cutoff q and 
the normalisation of the CR distribution function (Jubelgas et al. 
2006). Current cosmological radiative simulations that do not in- 
clude feedback from AGN bubbles (Sijacki & Springel 2006) have 
an over-cooling problem in the cluster centre where they overpro- 
duce the amount of stars. This problem is larger in cool core clus- 
ters which did not recently experience a major merger. Compared 
to X-ray observations, the simulated clusters show enhanced cen- 
tral gas densities, too small central temperatures, too strong central 
entropy plateaus, and an enhanced baryon fraction that is driven by 
the large stellar mass fraction towards the cluster centre. These are 
all manifestations of the cooling flow problem of simulated clus- 
ters, an issue that universally occurs in ACDM simulations, inde- 
pendent of the employed numerical method (e.g., Borgani et al. 
2004; Kravtsov et al. 2005; Ettori et al. 2006). 

In our plots of cluster profiles, the error bars represent the 
standard deviation <x from the sample mean of the indicated sub- 
sample of clusters. The error on the mean can easily be obtained 
by means of <rj VA' - 1, where A' is the number of clusters in the 
sample, i.e. N = 2 (small CC), N = 3 (large CC), N = 5 (all 
CC), N = 9 (merging cluster). The colour scheme for the profiles 
of our radiative simulations is black for our reference simulations 
without CRs, red for simulations where CRs are injected only at 
structure formation shocks using Mach number dependent scheme 
{CR-shock model), and blue for the complete model that accounts 
for CR injection from structure formation shocks on large scales 
(with the Mach number dependent scheme) and supernova shocks 
within the ISM of individual galaxies with a constant energy injec- 
tion efficiency of f S N = 0.3. Although our complete model might 
seem to be more relevant astrophysically, both CR models are cho- 
sen to represent extreme cases that are bracketing the realistic sce- 
nario in between for the following reasons. The CR-shock model 
approaches the realistic case from below because it does not take 
into account CRs from supernova shock waves. Our complete CR 
model brackets the model from above because it assumes a high CR 
injection efficiency from supernovae in addition to the too large star 
formation rate of current radiative simulations. This possibly over- 
estimates the CR energy injection into the ICM in the complete 
model. After describing the simulated profiles, we synthesise our 
results in Sect. 3.3.7 and provide further support for our interpreta- 
tion by considering CR related physical processes in more detail in 
Sect. 3.4. 



3.3.1 Density and temperature profiles 

Results for the appropriately scaled average density and tempera- 
ture profiles of our radiative simulations are shown in Fig. 5. Com- 
paring both CR simulations to our reference simulation, the most 
striking feature is the density enhancement in the central regions at 
r < 0.04i? vir . Our CR-shock model shows a density enhancement 
at larger radii which is strongest in small cool core clusters and vis- 
ible up to radii as large as r = 0.5 i? v ir- In contrast, our complete 
model shows a density decrement at intermediate radii that almost 
disappears in small cool core clusters. 

Similarly, the temperature in the central regions is reduced in 
our CR simulations while it almost coincides at larger radii with 
our reference simulation. The cluster-to-cluster variance is largest 
in our merger simulation due to ongoing virialisation processes me- 
diated by structure formation shock waves. 



3.3.2 Pressure profiles 

Figures 6 and 7 compare profiles of our sample of all cool core 
clusters (CC clusters, on the left-hand side) to our sample of merg- 
ing clusters (on the right-hand side). Our samples of small and large 
CC clusters differ only marginally for these profiles. While the ther- 
mal pressure profile of the CC clusters follows a cored /3-profile, 
the pressure profile of merging clusters almost shows a flat single 
power-law behaviour that is due to the ongoing merging activities 
(cf. Fig. 6, top panels). The thermal pressure of the CC clusters 
in our CR-shock model exceeds the one in our reference simula- 
tion, while the thermal pressure in the merging clusters shows no 
significant difference among our different models apart from the in- 
nermost central region where we observe an excess in both of our 
CR models. 

In contrast, the CR pressure profiles do not show a core; their 
profile increases inwards roughly as F C r k r' 1 up to the exter- 
nal regions of the ISM of the cD galaxy at r cc 0.04 R vil , where 
their profile steepens considerably and scales as P C r k r 2 m our 
CC clusters. The error bars of the CR pressure profiles have been 
suppressed for clarity. The CR pressure profiles in both of our CR 
models are similar to each other: the profile in the CR-shock model 
shows a smoother behaviour while the one in our complete model 
has a larger radial variation and a higher pressure level due to the 
additionally injected CR pressure from supernovae within galaxies. 
The hypothesis that the supernova CRs are irrelevant for the total 
pressure budget is however inadmissible because simulations that 
consider only CRs from supernovae shocks show very similar pres- 
sure profiles (Jubelgas et al. 2006). The reason lies rather in the fact 
that any CR population is processed by adiabatic loss and transport 
processes in a comparable fashion. 



3.3.3 Relative CR pressure profiles 

It is very instructive to study the relative CR pressure profiles 
XcrO") = Pcair) I P t h(r) in Fig. 6, middle panels. Generally, X C R(r) 
is larger in merging clusters compared to CC clusters and does 
not adopt values below X C r = 0.1 at any radius for the complete 
CR model in merging clusters. As already mentioned in Sect. 3.1, 
strong merger shock waves efficiently inject CRs at the exter- 
nal cluster regions and mix the highly CR-enriched intergalactic 
medium outside clusters with the ICM, thus boosting the relative 
CR pressure during a merger. This also results in large cluster-to- 
cluster variations at some given radius. Note that the error bars of 
the relative CR pressure profiles have been suppressed for clarity. 

Compared to the very smooth, dynamically unimportant pro- 
file of our non-radiative simulations, X CR (r) shows a steep increase 
towards the cluster centre within the ISM of the cD galaxy and 
much larger radial variation in these radiative simulations. This is 
due to the short thermal cooling time scales within dense galactic 
environments compared to the comparatively long CR cooling time 
scales (cf. Fig. 4). There, the thermal gas cools on a small timescale 
of t — 10 7 yr, diminishes pressure support, condenses out of the 
gaseous phase and eventually forms stars. In the process, the CR 
pressure reaches values that exceed the equipartition with the ther- 
mal pressure. In the outer parts of the cluster, X C R.(r) is a rising 
function with increasing radius due to the same reasons mentioned 
above in our non-radiative simulations. 
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Figure 5. Average profiles of our sample of clusters in our radiative simulations at redshift z = 0. Shown are the appropriately scaled profiles of the gas density 
p gas (left-hand side) and temperature T (right-hand side). Colour scheme: black: reference simulations without CRs, red: CRs are injected only through shock 
acceleration using Mach number dependent scheme, blue: complete CR-model including CR injection from supernovae. 



3.3.4 Effective adiabatic index 

In our CC cluster sample, the effective adiabatic index y cff departs 
from its canonical value y,h = 5/3 (in the absence of CRs) strongly 
in the cluster centre inside the ISM of the cD galaxy and at the 
external regions of the cluster where the relative CR pressure profile 
is dynamically important (cf. Fig. 6, bottom panels). The composite 
of CRs and thermal gas acquires there a higher compressibility that 
requires a larger overdensity in order to exert the same pressure 
that balances the gravitational pressure from the combination of 
dark matter and gas. The value of y c ff in our sample of merging 
clusters is everywhere softer compared to the canonical adiabatic 
index although there is a large variation among different clusters. 
This again reflects the higher CR pressure level in merging clusters 
due to the high CR injection efficiency of merger shock waves. 

Interestingly, the adiabatic index of the CR population ycR m 
our radiative simulations is harder compared to our non-radiative 
simulations. This can be traced back to the smaller overdensity in 
radiative simulations and thus the smaller CR cooling rate due to 
Coulomb interactions that re-thermalize the low-energy part of the 
CR momentum distribution function and eventually yield an ultra- 
relativistic CR population. The lower overdensity in our radiative 
simulations basically stems from the fact that gas mass from inter- 
mediate cluster scales is removed since it replenished the condens- 
ing star forming material inside galaxies and most importantly of 
the central cD galaxy. 

3.3.5 Baryon fraction 

Figure 7 shows the baryon fraction scaled with the universal value, 
b = /bar y (< x) Q. m /Q.]„ that reaches values of unity at the virial ra- 
dius independent of the adopted model and dynamical state of the 



cluster. The profile of b rises steeply towards the centre, which is 
mostly due to the increase of the stellar mass fraction relative to the 
dark matter. Additionally, Fig. 7 shows the gas fraction scaled with 
the universal value, / gas (< x)Q. m /Qb. Its profile is decreasing in- 
wards in contrast to the scaled baryon fraction. This shows clearly 
that the stellar mass fraction as well as the stellar mass-to-gas frac- 
tion is not constant as a function of radius but strongly increasing 
inwards. As a word of caution, a fraction of the inward increase of 
b is caused by the over-cooling problem in the cluster centre where 
we form too many stars that are too blue compared to observations. 
We note that this is a result of the well-known cooling flow prob- 
lem in simulated clusters, which is independent of the employed 
numerical method (e.g., Borgani et al. 2004; Kravtsov et al. 2005; 
Ettori et al. 2006, and references therein). CR physics with CR 
source terms studied in this paper such as structure formation shock 
waves and supernovae shocks on galactic scales seems to reinforce 
this problem instead of solving it. Eventually, self-regulating heat- 
ing sources related to AGNs at cluster centres might be required to 
solve this problem. 

The value of b in our sample of merging clusters is increased 
by roughly 10 per cent in our CR models compared to our reference 
simulation while the enhancement is smaller in the CC cluster sam- 
ple. The simulations show a decreasing gas mass fraction towards 
the centre for reasons that have been explained in Sect. 3.2.3. 



3.3.6 Cumulative stellar mass 

The profile of the cumulative stellar mass normalised by the virial 
mass of the cluster is shown in Fig. 7. One notes that there is a 
higher fraction of stellar mass formed relative to M y - n in our small 
CC clusters compared to the larger ones. This might be partly due 
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Figure 6. Average profiles of our sample of all cool core clusters (left-hand side) and merging clusters (right-hand side) in our radiative simulations at redshift 
Z = 0. Shown are the appropriately scaled profiles of the thermal and CR pressure (P t (, and Pqr, the latter with symbols), relative CR pressure Xcr = PcR/Pth* 
and the effective adiabatic indices y c ff (solid) in comparison with the CR adiabatic indices ycR (dashed). Colour scheme: black: reference simulations without 
CRs, red: CRs are injected only through shock acceleration using Mach number dependent scheme, blue: complete CR-model including CR injection from 
supernovae. 
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Figure 7. Average profiles of our sample of all cool core clusters (left-hand side) and merging clusters (right-hand side) in our radiative simulations at redshift 
7, = 0. Shown are the appropriately scaled profiles of the baryon fraction b (top panels) and the stellar mass M stars (< x) (bottom panels). Colour scheme: black: 
reference simulations without CRs, red: CRs are injected only through shock acceleration using Mach number dependent scheme, blue: complete CR-model 
including CR injection from supemovae. 



to the above mentioned over-cooling problem which has a higher 
relative impact in smaller systems. There are more stars formed 
in our complete CR model on scales r < 0.3W v ir irrespective of the 
dynamical state of the cluster while our CR-shock model shows less 
formed stellar mass in the outer cluster regions in our CC clusters. 
This shows again that the CR physics studied here aggravates the 
cooling flow problem in contemporary simulations. 



3.3.7 Discussion of CRs profiles in radiative simulations 

The CR distribution in the dilute outer cluster regions is dominated 
by the CR injection characteristics with its strong structure forma- 
tion shocks that efficiently accelerate CRs into the IGM. The rel- 
ative CR pressure X C r decreases as we move inwards due to the 
same reasons that apply in the case of non-radiative simulations: 
(1) weak central flow shocks are inefficient in accelerating CRs, (2) 
adiabatic compression of a composite of CRs and thermal gas dis- 
favours the CR pressure relative to the thermal pressure due to the 



softer equation of state of CRs, and (3) CR loss processes are more 
important at the dense centres. 

Rather than becoming dynamically unimportant as it is the 
case in our non-radiative simulations, the relative CR pressure in- 
creases strongly towards dense substructures if we consider radia- 
tive losses of the thermal gas due to the long CR cooling time scales 
compared to those of the thermal gas. The thermal gas cools much 
faster radiatively and diminishes its pressure support while it con- 
denses out and forms stars. The high fraction of pressure support 
provided by CRs yields a higher compressibility of the composite 
fluid that requires a larger overdensity in order to balance the grav- 
itational pressure from the dark matter and gas. This increases the 
density of each galactic substructure as well as in the cluster cen- 
tre and leads to a pressure enhancement, provided the hydrostatic 
equation is applicable, dP/(dr) = p gas GM(< r)r~ 2 . This higher 
density leads to a higher rate of star formation and thus increases 
the central baryon fraction since gas from larger scales has to re- 
plenish the condensed gas in order to maintain hydrostatic equilib- 
rium. 
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In the case of ongoing merger activity, the relative CR pressure 
is boosted due to a combination of two effects: (1) a merger event 
triggers large random motions and merger shock waves that steepen 
considerably while they are running into the shallower cluster po- 
tential where they efficiently inject CRs. (2) These large random 
motions mix the highly CR-enriched intergalactic medium outside 
clusters with the ICM. Weak virialisation shocks traversing the 
cluster after the merger thermalize random gas motions, thus in- 
crease the thermal energy and yield a decrease of X C r. CR loss pro- 
cesses and adiabatic compression furthermore decreases X C r in in- 
termediate cluster regions with a comparatively long thermal cool- 
ing time scale. 

3.4 Cosmic ray related physical processes 

In the following sections we study the physical processes related 
to CRs in more detail, which complements the previous sections. 
In particular, we will consider projected maps of physical cluster 
properties, and analyse the Mach number statistics during clus- 
ter formation. Finally, Sect. 3.5 studies the consequences of CRs 
for thermal observables such as X-ray emission and the Sunyaev- 
Zel'dovich effect. 

3.4.1 Maps of physical quantities 

We produced projected maps of the density, the mass-weighted 
temperature, CR related quantities, and thermal cluster observables. 
Generally, a three-dimensional scalar field a(r) along any ray was 
calculated by distributing the product of a(r) and the specific vol- 
ume M„lp a of the gas particles over a grid comoving with the cos- 
mic expansion. This yields the projected quantity A(r ± ): 

A(r ±Jj ) = -i-V%- W a jj(r ±J j - r a ), (4) 

where W a jj is the value of the projected smoothing kernel (nor- 
malised to unity for the pixels covered) of an SPH particle a at 
comoving grid position r ±> y, and L 2 ^ is the comoving area of the 
pixel. 

The surface mass density map was produced by projecting 
the quantity a a = p„/(ilbPcr) and the mass-weighted temperature 
map by setting a a = T a p a divided by the mass projection. Fig- 
ure 8 shows the radiative simulation of a super-cluster region (with 
the identifier g8) which is dominated by a large cool core clus- 
ter and surrounded by smaller satellite clusters. While the ICM 
of the central massive cluster reaches a temperature of the or- 
der of its virial temperature of kT vir = 13.1 keV, the surrounding 
warm-hot intergalactic medium (WHIM) acquires temperatures of 
kT ~ (10~ 2 - 10 _1 ) keV. The spatial distribution of shock strengths 
can be studied best by looking at the Mach numbers weighted by 
the energy dissipation rate at shocks. The colour hue of the lower 
panel on the left-hand side of Fig. 8 has been chosen to represent 
the projected quantity a„ = M a s^ s a divided by the projected en- 
ergy dissipation rate (£diss>ios, where M a denotes the Mach number 
of the SPH particle. The brightness scales logarithmically with the 
projected dissipation rate (fidisshos- Within this super-cluster region 
most of the energy is dissipated in weak internal shocks with Mach 
numbers M < 2 which are predominantly central flow shocks or 
merger shock waves traversing the cluster centre. Collapsed cos- 
mological structures are surrounded by several shells of external 
shocks with successively higher Mach numbers, but they play only 
a minor role in the energy balance of thermalization as can be in- 
ferred by its dim brightness. Clearly visible are spherical shells of 



shocks at different radii from the cluster centre. Two distinct accre- 
tion shocks at distances of 2 and 3 h~ [ Mpc to the cluster centre are 
visible, followed by shells of stronger shocks further outwards. 

Finally, the mass-weighted CR pressure relative to the total 
pressure X C r = Pcr/Pm was obtained by projecting a a = p a X C R, a 
divided by the mass projection (lower panel on the right-hand side 
of Fig. 8). The projections of the relative CR pressure support the 
picture inferred from our profiles. While X C r acquires compara- 
tively high values within the WHIM that are hydrodynamically 
important, their importance diminishes (on average) within each 
galaxy cluster for reasons laid out above (cf. Sect. 3.3.7). Within 
each individual galaxy, the CR pressure reaches equipartition or 
dominates the thermal pressure as can be seen by the numerous yel- 
low points sprinkled over the map, each corresponding to a galaxy. 
The time evolution of the relative CR pressure from z = 1 until 
today and the comparison of X C r for a large and a small CC clus- 
ter can be studied in more detail in Fig. 9, that shows the central 
regions of the galaxy cluster 1 and 10. At z = 0, the iso-contour 
of Xqr = 0.1 (the transition from green to red for our colour 
scale) roughly coincides with the virial radius in our small cluster 
R vil - 0.7 /r 1 Mpc while it approximately lies at a scale of twice 
the virial radius in the large cluster (R va ^ 2 /T 1 Mpc). This already 
hints that weaker shocks are responsible for thermal dissipation in 
high mass halos compared to low mass halos. We will study this 
effect in more detail in the next section. The second striking obser- 
vation is that at high redshifts the warmer colours dominate over 
the cooler ones, indicating a higher level of relative CR pressure 
at high redshift. This is because in the hierarchical model of struc- 
ture formation, the characteristic halo mass and thus the associ- 
ated characteristic temperature grow with cosmic time. The higher 
sound speeds of the diffuse gas at later times decreases the charac- 
teristic Mach numbers of the shocks responsible for thermalization 
and thus decreases the acceleration efficiency of CRs. 

3.4.2 Mach number statistics 

Earlier work on this subject (Miniati et al. 2000; Ryu et al. 2003; 
Pfrommer et al. 2006) classified structure formation shocks into 
two broad populations which are labelled as internal and external 
shocks. These categories have been distinguished on whether or not 
the associated pre-shock gas was previously shocked in the cosmic 
history. Since we are also interested in assessing the properties of 
shock waves in radiative simulations, it turns out to be favourable 
not to use a thermodynamical criterion such as the temperature but 
a criterion based on the overdensity 5 %m = p eas /(^bPcr a 3 ) - 1 in 
order not to confuse the shock definition in the presence of radia- 
tively cooling gas in galaxies (here, p gas denotes the physical gas 
mass density and a is the cosmic scale factor). In practice, we use 
the criterion of a critical pre-shock overdensity <5 gas > 10 or the 
physical gas density p eas > 7.5 x 10~ 24 g cm" 3 for the classification 
of an internal shock. With this criterion, external shock surfaces 
surround filaments, sheets, and haloes where the pristine adiabati- 
cally cooling gas is shocked for the first time while internal shocks 
are located within the regions bound by external shocks and are 
created by flow motions accompanying hierarchical structure for- 
mation. Their population includes accretion shocks produced by 
infalling material along the filaments into clusters, merger shocks 
resulting from infalling haloes, and flow shocks inside nonlinear 
structures which are excited by supersonic motions of subclumps. 
The thermal gas of the warm-hot intergalactic medium has a low 
sound velocity c = ^fy^Pjp = \jjihkT l(]j.m p ) owing to the low 
temperature, so once the diffuse gas breaks on mildly nonlinear 
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Figure 8. Environment of the largest galaxy cluster in our sample (a cool core cluster) in the radiative simulation including CRs from structure formation 
shocks only. Shown are the projected density, mass-weighted temperature, the Mach number of shocks weighted by the energy dissipation rate in colour (while 
the brightness displays the logarithm of the dissipation rate), and the mass weighted CR pressure relative to the total pressure. We performed line-of-sight 
averages to obtain the projected quantities. 



structures, strong shock waves develop that are characterised by 
high Mach numbers M = v s /c. Nevertheless, the energy dissipa- 
tion rate of internal shocks is always higher compared to external 
shocks within the WHIM because the mean shock speed and pre- 
shock gas densities are significantly larger for internal shocks. 

To quantify these considerations, we compute the differential 
Mach number distribution weighted by the shock-dissipated en- 
ergy d 2 £diss(tf, Af)/(dloga dlog M) at different redshifts charac- 
terised by scale factor a. The top left-hand panel of Fig. 10 shows 
this Mach number distribution in our non-radiative simulation of 
a super-cluster region (g8), while the top right-hand panel shows 
this distribution for the radiative simulation using the same initial 
conditions. The lower left-hand panel shows both distributions in- 
tegrated over the scale factor, d£ diss (A1)/(dlog M). Internal shocks 
are shown with dotted lines and external shocks with dashed lines. 



The lower right-hand panel shows again the integrated Mach num- 
ber distribution. This time, however, the shock population is sub- 
divided into two categories that are separated by the overdensity 
threshold <5 gas > 10 4 to study the origin of the flat horizontal tail of 
the Mach number distribution extending to very high Mach num- 
bers in our simulation with radiative gas physics. 

Several important points are apparent. (1) The non-radiative 
simulation shows a smooth distribution that is skewed towards 
lower values of the Mach number. The mean of the Mach number 
distribution weighted by the dissipated energy decreases as cosmic 
time evolves, i.e. the average shock becomes weaker at later times. 
(2) In contrast, the radiative simulation shows a broad symmetric 
distribution in logAl at early times and develops later on a peak 
at low Mach numbers that resembles that in the non-radiative case 
followed by a significant disturbed tail towards higher Mach num- 
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Figure 9. Time evolution of the mass-weighted CR pressure relative to the total pressure in the radiative simulation including CRs from structure formation 
shocks. Upper panels show the clusters at z, = 1, the lower panels at z = 0. Left panels: small cool core cluster with a virial mass of8.8xl0 13 /!" 1 M o atz = 0. 
Right panels: large cool core cluster with a virial mass of 1.8 X 10 15 hT l M e . 



bers of M > 10. Due to this tail, the mean of the distribution is not 
a monotonic function of time any more and achieves high values of 
<Md£ diss /(dlogM»/<d£ dlss /(dlogM)> - 100 already at z = 10. 
This reflects the merging history of the clusters forming in that 
simulation in combination with the ongoing development of dense 
cooling regions that permit the formation of strong shocks during 
merger events. (3) The time integrated Mach number distribution 
weighted by the dissipated energy, d£(ji ss (A4)/(dlog M), peaks at 
Mach numbers M < 2. Compared to full-size cosmological sim- 
ulations, the peak position remains unchanged while the smooth 
tail in our non-radiative super-cluster simulation is shifted towards 
lower values (cf. Pfrommer et al. 2006, Fig. 6). In non-radiative 
simulations, the peak of the Mach number distribution originates 
from weak flow shocks internal to dense halo regions, while the 
smooth tail towards strong shocks results from external shocks 



where the diffuse gas breaks on mildly nonlinear structures. These 
under-dense regions are under-represented in our super-cluster sim- 
ulation and thus we miss a subsample of this shock population com- 
pared to the cosmic mean. (4) The horizontal tail extending towards 
high Mach numbers in radiative simulations can unambiguously be 
attributed to internal shocks within cooled high-density regions of 
£ g as > 10 4 . Part of this result is due to the over-cooling of these 
dense regions in our radiative simulation, which yields an over- 
estimate of the Mach numbers due to the too small sound speeds. 
However, the main contribution in terms of energy dissipation orig- 
inates from internal shocks that occur in regions with a pre-shock 
overdensity 10 < <5 gas < 10 4 . This is because most of the ICM that 
is being virialised lies simply at overdensities S gas < 10 4 . (5) In both 
types of simulations, there is an increasing amount of energy dissi- 
pated at shock waves per logarithmic interval in the scale factor as 
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Figure 10. Influence of radiative cooling and star formation on the Mach number statistics of high-resolution cluster simulations (g8). The top left-hand panel 
shows the differential Mach number distribution d?E^(a, A1)/(dloga dlog M) for our non-radiative simulation while the top right-hand panel shows this 
distribution for the radiative simulation. The lower left-hand panel shows both distributions integrated over the scale factor, d£ ( ji ss (A1)/(dlog M). External 
shocks (occurring in overdensities <S gas = p gas /(flij p„ cT 3 ) — 1 < 10) are shown with dashed lines and internal shocks (<S gas > 10) with dotted lines. The lower 
right-hand panel shows again the integrated Mach number distribution, dEdi ss (A1)/(dlog M). Here, the shock populations are subdivided into two categories 
that are separated by the overdensity threshold (5 gas = 10 4 . This shows unambiguously that the flat tail of the Mach number distribution in radiative simulations 
originates from internal shocks within cooled high-density regions that are almost absent in non-radiative simulations. 



the universe evolves because the mean shock speed is significantly 
growing when the characteristic mass becomes larger with time. 
This trend starts to level off at redshift z—\ when galaxy clusters 
are in place. 

Figure 11 compares the shock strengths of our super-cluster 
region (g8) with that of an isolated small cluster region (g676). 
While the dissipated energy integrated over the cosmic forma- 
tion history in our simulation of the isolated small cluster (g676) 
amounts to E^ s = 2.4 x 10 62 erg, in our large super-cluster simula- 
tion (g8) the amount of dissipated energy was two orders of mag- 



nitude larger, £dj ss = 2.4 x 10 64 erg. Normalising the Mach number 
distribution d£' diss (A4)/(dlog M) with E iiS!i , we observe in our ra- 
diative as well as in our non-radiative simulation that the character- 
istic Mach number which is responsible for the energy dissipation 
decreases with cluster mass due to the increasing sound speed in 
more massive systems. In other words, weaker shocks are respon- 
sible for the energy dissipation in larger systems. Apparently, the 
characteristic shock speed in these large systems does not increase 
accordingly to make up for the increase in sound speed. This con- 
firms the suggestion of Jubelgas et al. (2006) that virialisation pro- 
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Figure 11. Influence of the cluster size on the Mach number statistics of high-resolution cluster simulations. The left-hand panel compares the appropriately 
scaled Mach number distribution d£ji ss (A1)/(dlogA1) of a large super-cluster region (g8) and a small isolated cluster region (g676) in our non-radiative 
simulation, while the right-hand panel shows the comparison in our simulation with radiative gas physics. Internal shocks (occurring at overdensities 5 ga s = 
/3gas/(£2b/°cr<z~ 3 ) > 10) are shown with dotted lines and external shocks with dashed lines. The total dissipated energy integrated over the cosmic formation 
history in our simulation of our isolated small cluster (g676) amounts to En ss = 2.4 X 10 62 erg, in our large super-cluster simulation (g8) the amount of 
dissipated energy was two orders of magnitude larger, E<ji ss = 2.4 X 10 64 erg. The characteristic Mach number which is responsible for the energy dissipation 
decreases with cluster mass due to the increasing sound speed in those systems. 
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Figure 12. Phase space distribution of the CR pressure relative to the thermal pressure at z = 0, Left-hand side: radiative simulation, right-hand side: 
non-radiative simulation of our small cluster 10. The black line shows a track of constant entropy of a composite of CRs and thermal gas during adiabatic 
compression (using the same normalisation in both panels to guide the eye). Dashed white line: star formation threshold of the simulation. Shown is the 
distribution of SPH particles with r < 10i? v j, that are identified by being in regions not contaminated by boundary dark matter particles in the simulation of 
cluster 10. 
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Figure 13. Average profiles of the Schwarzschild criterion for convective instability of all clusters. Left-hand side: non-radiative simulations, right-hand 
side: radiative simulations. Positive values of this criterion indicate convective stability. Colour scheme: black: reference simulations without CRs, red: CR are 
injected only through shock acceleration using Mach number dependent scheme, blue radiative: full model including CRs from supernovae, blue non-radiative: 
CR shock acceleration with constant energy injection efficiency. 



cesses in more massive systems is mediated by shocks with smaller 
Mach numbers. Since strong shocks are more effective in acceler- 
ating CRs, the relative CR pressure level is increased in low mass 
halos as shown in Table 3 and Fig. 9. 



3.4.3 Cosmic ray phase space distributions 

Figure 12 compares the phase space distribution of the CR popu- 
lation at i = of a radiative simulation to a non-radiative simula- 
tion. Specifically, we present the CR pressure relative to the ther- 
mal pressure, X C r, versus the overdensity of the gas (5 aas + 1 = 
Pga S /(fJbPcr<3 3 )- Thus, this figure shows the variance in Xqr at 
a given overdensity 6 giis in our simulations of cluster 10. At low 
densities, both phase space distributions resemble each other. Low 
density regions are characterised by strong shocks that efficiently 
accelerate CRs and yield a high value of Xqr. As we move in- 
wards the cluster, the density increases and the characteristic shock 
strength decreases which makes CR acceleration less efficient; an 
effect that reduces Xcr. Adiabatic compression of a composite of 
CRs and thermal gas disfavours the CR pressure relative to the ther- 
mal pressure due to the softer equation of state of CRs. This is vi- 
sualised by the black line which shows a track of constant entropy 
X C r oc (5gaf 7lh of a composite of CRs and thermal gas during adi- 
abatic compression. Finally, CR loss processes are more important 
in denser regions, a third process that reduces Xcr. While the story 
ends here for the non-radiative simulations, it goes an important 
step beyond in radiative simulations. There, the thermal gas cools 
on a small timescale, diminishes pressure support, condenses out 
of the gaseous phase and eventually forms stars at the characteris- 
tic density threshold as indicated by the white dashed line. Mean- 
while, the relative CR pressure increases due to the comparatively 
long CR cooling time leading to a CR pressure that dominates over 
the thermal one at densities S glis > 10 4 . Finally, at the star formation 
threshold, we obtain a very broad distribution of Xcr. 



3.4.4 CR buoyancy - the generalised Schwarzschild criterion 

Figure 13 shows the generalised Schwarzschild criterion for con- 
vective instability of our sample of all clusters. The criterion for 
convective stability is given by 



A In A 
Ar 



> 0, 



(5) 



where the effective entropic function is denoted A = (P th + 
Pcr) p~ 7M (cf. Eqn. (A8)). The ICM of our non-radiative and radia- 
tive simulations is overall convectively stable. This is no surprise 
since in any convective unstable region, the dynamics of the insta- 
bility will quickly restore a stable configuration, thus instabilities 
are self-erasing. The radiative simulations show a larger cluster- 
to-cluster variation with respect to individual radial bins, prefer- 
ably in merging clusters that are formally convectively unstable. 
Infalling clumps of matter distort the spherically averaged profile 
such that the Schwarzschild criterion classifies them as convec- 
tively unstable. The simulations including CR physics do not show 
significantly different behaviour compared to our reference simula- 
tions. CRs from supernovae shocks that are injected near the cluster 
centre could have triggered convective instabilities as proposed by 
Chandran (2004), but this is not really seen in our simulations. The 
situation might be different and more unstable if one included an 
additional central source in the form of CRs that are diffusing out 
of AGN bubbles at cluster centres. 



3.5 Thermal cluster observables 

3.5.1 Enhanced X-ray emission 

Our map of the bolometric X-ray surface brightness in the 
bremsstrahlung regime was produced by means of Eqn. (4), pro- 
jecting the quantity 



Ao 



kT \ 1/2 / o x2 



keV 



pm p 



(6) 
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Small cool core cluster: Large cool core cluster: 
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Figure 14. Difference of the X-ray surface brightness in a radiative simulation with CRs and the corresponding reference simulation without CRs. Left-hand 
side: small cool core cluster with a virial mass of 8.8 X 10 Mpc at z = 0, right-hand side: large cool core cluster with a virial mass of 1.8 X 10 15 h Mpc. 
Upper panels: simulations only account for CR acceleration at structure formation shocks, lower panels: CRs from supernovae are also considered. The relative 
difference of the integrated X-ray luminosity (excluding the central over-cooling region of r < 5 h~ 1 kpc) is given in the inlay. 



/27rkeV\ 1/2 2 s ne 6 , 1 +X S 

Ao = hi Ti -SB(T a )^ 2 — ^ (7) 

\ 3 m c / 3/1 m s c- 1 

= 1.78 x 1(T 24 erg cm" 3 s" 1 . (8) 

Here, T a denotes the temperature of the SPH particle (labelled 
with a), h denotes the Planck constant, e the elementary charge, 
X H = 0.76 is the primordial hydrogen mass fraction, /i = 4/(5X H + 
3) = 0.588 is the mean molecular weight of fully ionised gas, and 
gs — 1-2 is the frequency- and velocity-averaged Gaunt factor (Ry- 
bicki & Lightman 1979; Spitzer 1978). To be self-consistent within 
the bremsstrahlung assumption, we assume a fully ionised gas of 
primordial element composition which consists of hydrogen and 
helium only. 

Figure 14 shows the difference of the X-ray surface brightness 



in a radiative simulation with CRs and the corresponding reference 
simulation without CRs. The small blue and yellow circles repre- 
sent individual galaxies in each of the simulations whose orbits did 
not match up when one includes different physical CR processes 
such that small changes in the energy budget modify the galactic 
orbits, accumulate and eventually lead to macroscopic separation 
of equivalent galaxies in the final spatial distribution. The left-hand 
side shows a small cool core cluster with a virial mass today of 
8.8 x 10 13 hT l Mpc (cluster 10) and the right-hand side a large cool 
core cluster with a virial mass today of 1.8 x 10 /T 1 Mpc (clus- 
ter 1). While the upper panels show the simulation that only ac- 
counts for CR acceleration at structure formation shocks, the simu- 
lations displayed in the lower panels take also CRs from supernovae 
into account. The relative difference of the integrated X-ray lumi- 
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Table 4. : CR induced relative differences of the X-ray luminosity and the integrated Compton->' parameter in a representative cross section of our cluster 

SAMPLE. 



CI. 


sim.'s 


dyn. state' 1 ' 


T m 

200 
[keV] 


AL X /Lf 
shock-CRs (3) 


AL x /lf 
compl. CRs (3) 


AY/Y 
shock-CRs (3) 


AY/Y 
compl. CRs (3) 


Ayo/yo 

shock-CRs <3) 


A>'o/yo 
compl. CRs (3) 


>'o 

[icr 5 ] 


1 


g8a 


CC 


13.1 


+23% 


+22% 


+0.1% 


-0.6% 


+ 8% 


±0% 


26.7 


2 


gla 


cc 


10.6 


+ 18% 


+ 14% 


+ 1.1% 


-0.1% 


+ 8% 


-4% 


15.4 


3 


g72a 
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9.4 


+ 9% 


+ 15% 


-1.6% 


-0.1% 


-7% 


-7% 


14.6 


5 


gib 


M 


4.7 


+ 16% 


+69% 


-1.5% 


+ 1.1% 


+ 7% 


-4% 


2.73 


9 


gld 


M 


1.7 


-10% 


+ 8% 


-4.5% 


-2.2% 


-12% 


- 1% 


0.475 


10 


g676 


CC 


1.7 


+29% 


+52% 


-0.8% 


-4.5% 


+23% 


+41% 


0.960 


11 


g914 


cc 


1.6 


+36% 


+23% 


-0.1% 


-3.7% 


+25% 


+ 5% 


1.01 



Notes: 

(1) The definitions for the dynamical state of the cluster and the virial temperature is given in Table 2. 

(2) The relative difference of the integrated X-ray luminosity has been computed without taking into account the central over-cooling region of r < 5 A -1 kpc. 
Positive values for the relative difference indicate an enhancement of the thermal cluster observable when taking CRs into account. 

(3) The shock-CR model is represented by our radiative simulations that account only for structure formation CRs while the complete model is based on our 
radiative simulations that take additionally supernova CRs into account. 



nosity (without the central over-cooling region of r < Sh' 1 kpc) 
is given in the inlay, appropriately colour-coded such that the in- 
tegrated CR effect on the X-ray luminosity can be easily read off 
from the colour bar. All these simulations show an enhancement 
of the central and integrated X-ray luminosity in our CR models 
compared to our reference models without CRs. This effect derives 
from the softer effective adiabatic index in the CR simulations due 
to the large ratio of the CR-to-thermal cooling time in the cluster 
centre, leading to a more compressible gas with an associated den- 
sity increase that is required to balance the gravitational pressure 
of the dark matter. The resulting enhancement of the X-ray sur- 
face brightness is more pronounced in small CC clusters due to the 
higher level of relative CR pressure in combination with a relaxed 
spherically symmetric structure. 

Our previous considerations are confirmed by analysing a rep- 
resentative subsample of our cluster simulations with respect to 
mass and dynamical state in Table 4. Positive values for the relative 
difference indicate an enhancement of the thermal cluster observ- 
able when taking CRs into account. While the differences of the 
cool core cluster simulations yield a consistent picture in our runs 
which account for CRs from structure formation shocks only, the 
sign and the magnitude of the CR-induced effects in merging sys- 
tems do not follow a uniform trend due to the merger induced inho- 
mogeneities in the ICM. In our complete CR model, the cluster-to- 
cluster variation of the CR enhancement of equally sized clusters is 
larger due the CR injection by supernovae, which traces the current 
star formation rate rather than being an integrated measure over the 
past cluster/core activity like in our CR-shock runs. 

We can get a quantitative idea of the approximate size of this 
CR effect by employing a toy model of cluster formation in the 
presence of a CR component which we compare to a model with- 
out CRs. Consider at some early stage around z ~ 2- 3 a composite 
of thermal gas and CRs that has been shock heated for the first time 
while it was collapsing into mildly non-linear structures. Today's 
analogue for such a state of the gas might be given by the WHIM 
which is characterised by a comparatively high relative CR pres- 
sure contribution of roughly X C r ~ 0.3 - 0.5. Let us assume that 
the initial thermal pressure in our reference model without CRs is 
equivalent to the initial total pressure in our composite model with 
CRs, P' lh = P thp + -P C r,o = P<h,o (1 + ^cr,o), where our reference 
quantities are denoted by a star. Neglecting radiative cooling and 
CR diffusion for simplicity, further accretion of new shells of mass 



onto the dark matter halo exerts a pressure that adiabatically com- 
presses the gas internal to the halo, yielding 

Ah = Ah,oC™ and P cr = Pcr,oC cr . (9) 

While the gas in our reference model is compressed by the factor 
C*, the compression factor in the CR model is given by C. Re- 
quiring that the final pressure after adiabatic compression in both 
models is equivalent to each other, we obtain the condition 

/ C \ y ' h l C \ rcR 

(75) + *cr,o(^) C«™-(1+Xciy)) = 0. (10) 

This equation can be solved for the combination C/C*. Assuming 
an initial relative CR pressure of Xcr,o — 0.3, a compression factor 
in the reference model of C* - 10 3 , and equivalent initial densities 
in both models, we obtain a density enhancement for our toy cluster 
in our CR model of Ap/p* =s 15 per cent relative to the reference 
model. The associated X-ray emission would be enhanced in the 
presence of CRs by a factor of AL X /L* X ^ 32 per cent compared 
to the model without CRs. This result should be taken with a grain 
of salt since we did not consider shock heating, radiative cooling, 
nor any CR transport and loss processes. It merely illustrates the 
important effect of the softer equation of state of CRs. 



3. 5. 2 Modified Sunyaev-Zel 'dovich effect 

The amplitude of the thermal Sunyaev-Zel'dovich effect is given 
by the Compton-y parameter and is proportional to the line-of-sight 
integral of the thermal pressure which is obtained by replacing a in 
Eqn. (4) by 

cr T 

a„ = (y lh - 1) X H %z, a H u a p a , (11) 

m e c 2 

where Xyi = 0.76 is the primordial hydrogen mass fraction, x e , is 
the ratio of electron and hydrogen number densities of the SPH par- 
ticle (labelled with a) which we dynamically track in our radiative 
simulations, p = 4/(3X H + 1 + 4X H x eja ) denotes the mean molec- 
ular weight of partly ionised gas, u a is the internal energy per unit 
mass, and p„ is the physical mass density of the SPH particle. 

Figure 15 shows the difference between the Compton-y pa- 
rameter in a radiative simulation with CRs and the corresponding 
reference simulation without CRs. The left-hand side shows a small 
cool core cluster with a virial mass today of 8.8 x 10' 3 /T 1 Mpc 
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Figure 15. Difference of the Compton-y parameter in a radiative simulation with CRs and the corresponding reference simulation without CRs. Left-hand side: 
small cool core cluster with a virial mass of 8.8 X 10' 3 h Mpc at z = 0, right-hand side: large cool core cluster with a virial mass of 1.8 X 10 15 ft -1 Mpc. Upper 
panels: simulations only account for CR acceleration at structure formation shocks, lower panels: additionally, CRs are also injected by supernova shocks. The 
relative difference of the integrated Compton-y parameter is given in the inlay. 



(cluster 10) and the right-hand side displays a large cool core clus- 
ter with a virial mass today of 1.8 x 10 15 /T 1 Mpc (cluster 1). While 
the upper panels show the simulations that only account for CR ac- 
celeration at structure formation shocks, the simulations displayed 
in the lower panels take also CRs from supernovae into account. 
The relative difference of the integrated Compton-y parameter is 
given in the inlay, appropriately colour-coded such that the CR ef- 
fect on the integrated Compton-y parameter can easily be read off 
from the colour bar. In the small cluster simulation, CR feedback 
leads to an enhancement of the central Compton-y parameter fol- 
lowed by a ring of a (much weaker) decrement. The more massive 
system shows a more asymmetric difference map with a larger con- 
trast at the centre that is roughly one order of magnitude larger 
compared to the smaller system, as expected for the larger virial 



mass. As analysed in Sect. 3.3, the high fraction of pressure sup- 
port provided by CRs in the central regions results in a higher com- 
pressibility of the composite fluid, which in turn requires a larger 
overdensity in order to balance the gravitational pressure from the 
dark matter and gas. This increases the central density of the clus- 
ter and leads to a pressure enhancement, provided the hydrostatic 
equation is applicable, dP/(dr) = p gas GM(< r)r~ 2 . Hence, the 
Compton-y parameter in our CR simulations shows a profile that is 
more peaked compared to our reference case. The increased cen- 
tral Compton-y parameter comprises a larger area in our CR-shock 
model compared to the complete CR model. This is due to the addi- 
tional CR injection from SN in our complete CR model that reduces 
the amount of thermal SN feedback and thus reduces the Compton- 
y parameter. 
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Table 4 confirms this finding in a representative cross section 
of our clusters by showing the central Compton-y parameter y as 
well as the relative difference of y . In practice, this quantity will 
be smoothed by the experimental beam which we did not try to 
model here. The behaviour of Ay /yo is much more regular in our 
shock-CR model since the CR effects are integrated over a longer 
timescale. In contrast, there is a larger variation of Ay /y in our 
complete CR model, even among nearly similarly evolving clus- 
ters (cluster 10 and cluster 11), due to the central CR injection by 
supernovae which reflects the current central star formation rate. 

The total Compton-y parameter Y of a galaxy cluster is given 
by the Compton-y parameter integrated over the solid angle element 
subtended by the cluster, 

Y = D 1 mg f dny(0)= f d 2 *y(*J = (12) 
j£2 ctaMr «^W. Tth - i m e c 

where Z) ane is the angular diameter distance to the cluster. Interest- 
ingly, the total Compton-y parameter seems to be robust with re- 
spect to CR feedback because it reflects the cluster's total thermal 
energy E th . We confirm this finding by looking at a representative 
cross section of clusters in our sample in Table 4. In large CC clus- 
ters, CR feedback has a negligible effect on the total Compton-y 
while in merging clusters, the total Compton-y parameter can ei- 
ther be slightly enhanced or decreased, depending on the merger 
virialisation progress of the induced shock waves leading to tem- 
porary inhomogeneities in the ICM. In small CC clusters, the value 
of Y is only slightly decreased in our CR-shock model while this 
effect is more pronounced in our complete CR model where Y is 
systematically decreased by 4 per cent due to the additional CR 
injection from SN that reduces the Compton-y parameter. 



4 SUMMARY 

We performed high-resolution simulations of a sample of 14 
galaxy clusters that span a mass range from 5 x 10 13 /r' M Q to 
2 x 10 15 /r' M to study the effects of cosmic rays (CRs) on ther- 
mal cluster observables such as X-ray emission and the Sunyaev- 
Zel'dovich effect. In a systematic study, we analyse the CR effects 
on the intra-cluster medium while simultaneously taking into ac- 
count the cluster's dynamical state as well as the mass of the clus- 
ter. The modelling of the cosmic ray physics includes adiabatic 
CR transport processes such as compression and rarefaction, and a 
number of physical source and sink terms which modify the cosmic 
ray pressure of each particle. The most important sources consid- 
ered are injection by supernovae (in our radiative simulations) and 
diffusive shock acceleration at cosmological structure formation 
shocks, while the primary sinks are thermalization by Coulomb in- 
teractions, and catastrophic losses by hadronic interactions. While 
the relative pressure contained in CRs within the virial radius is of 
the order of 2 per cent averaged over our cluster sample in non- 
radiative simulations, their average contribution rises to 32 per cent 
in our simulations with dissipative gas physics including radiative 
cooling, star formation, and supernova feedback. 

Our main findings for the CR pressure relative to the thermal 
pressure X C r = Pcr/P± can be summarised as follows. The CR 
distribution in the dilute outskirts of cluster is dominated by the ef- 
fective CR acceleration at strong structure formation shocks such 
as accretion and merger shocks. The relative CR pressure Xqr de- 
creases as we move inwards to smaller cluster radii that are char- 
acterised by an increase of the sound speed and the density of the 
intracluster medium (ICM): (1) weak central flow shocks are ineffi- 



cient in accelerating CRs, (2) adiabatic compression of a composite 
of CRs and thermal gas disfavours the CR pressure relative to the 
thermal pressure due to the softer equation of state of CRs, and (3) 
CR loss processes are more important at the dense centres. 

Interestingly, if we include radiative losses of the thermal gas, 
the relative CR pressure increases strongly towards the cluster cen- 
tre and inside dense galactic substructures within the ICM, because 
of the long CR cooling times compared to those of the thermal gas 
in these cooling environments. The thermal gas cools much faster 
radiatively and diminishes its pressure support while it condenses 
out and forms stars that do not participate in the pressure balance. 
In contrast, at a fixed density, the CR cooling time of an aged CR 
population remains almost constant as the thermal gas is cooling. 
High fraction of pressure support provided by CRs yields a higher 
compressibility of the composite fluid that requires a larger over- 
density in order to balance the gravitational pressure from the dark 
matter and gas. This increases the density of each galactic substruc- 
ture as well as in the cluster centre where it leads to a pressure en- 
hancement, provided the hydrostatic equation is applicable as it is 
the case in cool core clusters. This higher density leads to a higher 
rate of star formation and thus increases the central baryon fraction 
strongly because gas from larger scales has to replenish the con- 
densed gas in order to maintain hydrostatic equilibrium. We note 
that in real galaxy clusters this effect is expected to be somewhat 
weaker, because here strong cooling flows are suppressed by some 
non-gravitational heating source. Nevertheless, this shows that CR 
populations injected from supernovae shocks on galactic scales and 
by structure formation shocks reinforce the cooling flow problem in 
contemporary simulations rather than solving it. CRs from super- 
novae shocks that are injected near the cluster centre could have in 
principle triggered convective instabilities in the lighter relativistic 
CR component, but this effect is not really seen in our simulations. 

The dynamical state of the cluster crucially determines the rel- 
ative CR pressure X C r'- (1) In the case of ongoing merger activ- 
ity, the relative CR pressure is boosted due to strong merger shock 
waves that effectively inject CRs and mix the highly CR-enriched 
intergalactic medium outside clusters with the ICM. (2) In a post- 
merger state, weak virialisation shocks traversing the cluster after 
the merger thermalize random gas motions, thereby increasing the 
thermal energy and yielding a decrease of X C r- CR loss processes 
and adiabatic compression furthermore decrease X C r in intermedi- 
ate cluster regions with a comparatively long thermal cooling time 
scale. (3) Small cool core clusters are characterised by a higher 
value of Xcr than more massive clusters, albeit the integrated Xqr 
is significantly reduced compared to our merging systems. This is 
mainly due to the nature of the virialisation process where weaker 
shocks are responsible for the energy dissipation in larger systems. 
This leads to a more efficient CR energy injection in small sys- 
tems and hence to a smaller value of X C r- Altogether, we observe a 
strong correlation between merging activity and the level of Xcr. 

How does a CR component affect thermal cluster observables? 
We found out that central cluster regions of low-mass clusters are 
most strongly affected by the CR component which leads to an 
enhanced central density and thermal pressure. Since the X-ray 
emission scales as the square of the gas density, the X-ray lu- 
minosity in low mass clusters of our CR simulations is boosted 
by up to AL X /L X 40 per cent. The amplitude of the thermal 
Sunyaev-Zel'dovich effect is also modified by CR feedback. While 
the central Compton-y parameter in these systems is increased by 
Ayo/yo - 25 per cent, the integrated Sunyaev-Zel'dovich effect ap- 
pears to be remarkably robust and the total flux decrement only 
slightly reduced by typically 2 per cent. In merging clusters, the 
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sign and the magnitude of the CR-induced effects are quite irregu- 
lar due to the merger induced inhomogeneities in the ICM. 
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APPENDIX A: SCHWARZSCHILD CRITERION FOR 
CONVECTIVE STABILITY 

This section follows the derivation by Chandran (2004) and derives 
a dimensionless criterion for convective stability in the presence 
of a composite fluid of thermal gas and CRs. Consider a parcel 
of thermal gas and CRs initially at a distance r = r from the 
cluster centre. The parcel is displaced by a radial distance Ar to 
r = r\ = ro + Ar. For simplicity thermal conduction and CR diffu- 
sion into and out of the parcel are neglected so that the parcel ex- 
pands adiabatically from the gas density plop'. At the final stage, 
the parcel is mixed with its surroundings. The volume occupied by 
the thermal gas and that occupied by the CRs expands by the same 
amount, implying 

p - = = i + S . (Al) 

P0 PCR.O 

The values of p, P th , and P C r in the parcel initially at the same av- 
erage values at r = ro which are denoted by po, _P t h,o, and Pcr.o- 
After the parcel is displaced radially outward, the new fluid quan- 
tities within the parcel are denoted p', P[^ and P'q^- The average 
fluid quantities (outside the parcel) at r = r\ are denoted p t , P^j, 
and Pcr.i- The total pressure is given by the sum of the partial 
pressures and denoted without any subscript, P = P th + P C r- The 
difference between the thermal gas density in the displaced parcel 
and its surroundings at r = r { is denoted by Ap = p' - p\. Other 
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fluid quantities such as the thermal pressure are denoted similarly, 
AP th = P' th - P th i, etc. Since the parcel expands adiabatically, we 
have 

P' ik = P< h , (l + 6y°> and P^ R = P C R,o(l+<5) rcR . (A2) 

Assuming that the turbulent velocities are subsonic, the total pres- 
sure inside the parcel approximately matches the average value out- 
side the parcel, P l = P'. We assume a small displacement |Ar| <k r 
and expand the average values of the density and the pressure out- 
side the parcel at r = n into their Taylor series, 

dp dP 
pi =s p + -f- Ar and P, =s P + — Ar. (A3) 
dr dr 

Employing adiabatic expansion of the total pressure within the par- 
cel (Eqn. A2), to lowest order in Ar/r we arrive at the following set 
of equations, 

1 dP 

— Ar, (A4) 



TrffP dr 

A P _ I 1 dP_ldp^ (A5) 



p \7tffP dr p dr 

AA _ (ldP_^dpU (A6) 



A \ P dr P dr) 

Here, we introduced for convenience the effective adiabatic index 
of the composite fluid, 

_ d logCPth + P C r) Til. Ah + 7cr -Pcr , A « 

dlogp Pth + PcR 

and the effective entropic function A = (P, h + PcR)p~ rdr - Since 
we kept only the lowest orders, it is not necessary to specify in 
Eqns. (A4) to (A6) whether the hydrodynamic quantities are eval- 
uated at ro or r\, and hence we dropped the subscripts on p, P, 
etc. Since the displacement Ar is a signed quantity, the criterion for 
convective stability is 

^ > 0. (A8) 
Ar 

For negligible CR pressure contribution, this criterion is equivalent 

to the classical Schwarzschild criterion 

AlnA th 

> 0. (A9) 

Ar 

Interestingly, the ICM can be convectively unstable even when the 
specific entropy of the thermal gas increases outward, as long as 
d_P C R/dr is negative and its absolute value is sufficiently large! 

This paper has been typeset from a TgX/ KTpX file prepared by the 
author. 
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